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ABSTRACT 



Context. Prominence seismology exploits our knowledge of the linear eigenoscillations for representative magnetohydro- 
dynamic models of filaments. To date, highly idealized models for prominences have been used, especially with respect 
to the overall magnetic configurations. 

Aims. We initiate a more systematic survey of filament wave modes, where we consider full multi-dimensional models 
with twisted magnetic fields representative of the surrounding magnetic flux rope. This requires the ability to compute 
accurate 2.5 dimensional magnetohydrodynamic equilibria that balance Lorentz forces, gravity, and pressure gradients, 
while containing density enhancements (static or in motion). 

Methods. The governing extended Grad-Shafranov equation is discussed, along with an analytic prediction for circular 
flux ropes for the Shafranov shift of the central magnetic axis due to gravity. Numerical equilibria are computed with 
a finite element-based code, demonstrating fourth order accuracy on an explicitly known, non-trivial test case. 
Results. The code is then used to construct more realistic prominence equilibria, for all three possible choices of a free 
flux-function. We quantify the influence of gravity, and generate cool condensations in hot cavities, as well as multi- 
layered prominences. 

Conclusions. The internal flux rope equilibria computed here have the prerequisite numerical accuracy to allow a yet more 
advanced analysis of the complete spectrum of linear magnetohydrodynamic perturbations, as will be demonstrated in 
the companion paper. 

Key words, solar physics, solar prominences - Instabilities - Magnetohydrodynamics (MHD) - Plasmas 
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1. Prominence seismology and equilibrium 
configurations 

One of the most fascinating phenomena embedded in the 
million degree solar coronal plasma is the presence of 
so-called filaments, which are plasma concentrations sus- 
pended magnetically against the downward pull of the so- 
lar gravitational field. They are up to 100 times colder and 
denser than their immediate surroundings, and in Ha obser- 
vations of the solar disk appear as dark features extending 
over huge distances. High resolution observations indicate 
that the filament is actually composed of many individual 
threads, down to the resolution l imit of modern observa- 
tions (about 100 km in width, see iLin et al.l (|2005T) ). while 
their length can reach several tens of megameters. When 
viewed at the solar limb, the filaments can be identified as 
prominences, and the surrounding magnetic geometry in- 
troduces a classification supporte d by early an alytical mag- 
netohydrostatic models (see e.g. lPriestill988l ): normal and 
inverse polarity prominences differ in the relative orienta- 
tion of the filament-carrying flux rope with respect to the 
underlying (and overarching arcade) magnetic orientation. 
Comparing their locations to photospheric magnetograms, 
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the filament threads appear primarily suspended above a 
magnetic neutral line. The study of solar filaments still 
poses many contemporary challenges to theoretical solar 
physicists, in terms of their sudden formation, their poten- 
tial to survive for months on end, and their small as well as 
large-scale dynamics. Reviews of the physics of solar promi- 
nences include one by lLabrosse et al.l (|2010l ). providing an 
overview of the prominence/filament thermodynamic pa- 
rameters o btained by modern sp ectroscopic inversions, and 
another by Mac kav et al.l (|2010t l. which compiled insights 
from essentially magnetohydrodynamic modeling that were 
compared with observations. 

A particularly active research area is known as promi- 
nence seismology, whose ultimate goal is the use of observed 
oscillatory motions in filaments to deduce their internal 
properties, by matching computed with detected frequen- 
cies, damping rates, and possibly spatio-temporal eigen- 
structure of t he vibration modes. As part of this effort, 
the review by iMackav et al.l (|2010fl collects insights from 
observations and theory, and points out how r ecent works 
accou nt for the multi-threaded fine structure (jLuna et al.l 
I2010D . as well as for non-adiabatic modifications to the wave 
m odes in prominences (ISoler et all 120081. Another review 
by lArregui fc BallesterJ ( 2011 ) concentrates on how the lin- 
ear MHD modes can help in explaining the observed at- 
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tenuation. To make analytic progress, virtually all of these 
efforts make dramatic simplifications about the overall ge- 
ometry of the magnetic field and thread characteristics: 
cylindrical flux tubes or segments embedded in uniform, 
thread-aligned magnetic fields form a recurring ingredient. 
When accounting for internal structure, density variations 
can explain the observed damping of global eigenoscillations 
in terms of resonant absorption and/or partial ionization, 
where global kink modes coup le to flux-surface localized 
Alfven modes (jSoler et al.ll2010D . Despite the model restric- 
tions, various indicative relations between thread-corona 
density contrast on the one hand, and mode damping rate 
to period ratio on the other hand, have been established 
by mean s of combined analytic al and numerical techniques 
(see e.g. iGoossens et ahl ([201 ll ) 1 ). At the same time, it is as 
yet unclear whether these relations generalize to the more 
realistic topologies where magnetic shear, gravity and pres- 
sure variations are all incorporated. 

A major obstacle on the road to generalizing these find- 
ings to more realistic prominence models, is that generating 
accurate equilibrium models for prominence-bearing flux 
ropes becomes a non-linear problem of its own. Even when 
we restrict the problem to a 2.5D magnetohydrostatic con- 
figuration by assuming translational invariance along the 
prominence axis, relatively scarce cases exist where an an- 
alytic sol ution is known . Notable exceptions form the early 
work by iDungevi (|1953ft , where for isothermal conditions 
the static force balance between the Lorentz force, the (uni- 
form) gravit y, and the pres s ure g radients could be solved 
analytically. lLow k Zhaiiel f|2004T ) developed under poly- 
tropic pressure-density relations, magnetostatic solutions 
for both normal and inverse polarity quiescent prominences. 
These latter authors pointed out how their construction 
avoids the free-boundary problem posed by matching the 
internal flux-rope solution, to the external (nearly poten- 
tial) atmospheric field. In practice, this was achieved by 
prescribing the prominence to be embedded in a horizon- 
tal, circular cylinder, whose boun dary then coincides with 
the internal-external transition. In lPetrie et al.l (|2007f) . this 
analytic approach was complemented and generalized with 
a numerical treatment where again both polarity types were 
realized, by iteratively solving the governing elliptic PDE 
which does fully account for the free-boundary problem. It 
was then demonstrated that a solution can be constructed 
with the characteristic, observed three-part structure of a 
cool dense prominence, surrounded by a cavity, within a 
flux rope in a hot coronal environment. 

We here revisit the problem of generating accurate, nu- 
merical equilibria that are representative of helical flux 
ropes, supporting dense prominences, possibly showing ax- 
ial (shear) flows. The need to demonstrate high accu- 
racy for the obtained numerical equilibria is a prerequi- 
site for any further analysis quantifying the magnetohy- 
drodynamic wave modes. This is common knowledge in 
fusion-rel ated plasma configurat ions, where MHD spec- 
troscopy (Gocdbloed et al.l Il993ft needs to take full ac- 
count of the detailed magnetic geometry to chart out, 
and eventually invert from, the complete MHD spectrum 
of waves and instabilities. Meanwhile, this approach has 
already bee n transported to astrophysically relevant ac - 
cretion tori (iGoedbloed et aU 12004 iBlokland et al.l 1 20071V 
where a toroidally caged, axisymmetric plasma torus has 
both toroidal (Keplerian-type) and possibly also poloidal 
flows in a delicate force balance with a central external 



gravitational field. The MHD spectroscopic determination 
of the eigenoscillations of these accretion tori brought out 
how flux-surface localized eigenmodes may become un- 
stable as a result of intricate mode couplings resulting 
from the poloidal cross-sectional variation. We here at- 
tempt to demonstrate the same ability to chart out all 
eigenmodes for realistically structured prominences, start- 
ing with the determination of the continuous parts of the 
MHD eigenspectrum. These MHD continua form the ba- 
sic organizin g structure of the entire c o llection of MHD 
eigen modes ([Goedbloed k Poedtsl 12004 iGoedbloed et al.l 
|2010D . and because of the multi-dimensional nature of the 
prominence-containing flux rope, will contain avoided c ross- 
ings and frequency gaps (|Poedts k Goossensl Il991h . To 
quantify these mode couplings, and to subsequently take 
the next step towards full diagnosis of all global discrete 
eigenmodes, we must ensure that our equilibrium solution 
is amenable to such an analysis. 

This paper is therefore organized as follows. We begin in 
Sect. [5] with reducing the governing equations for a transla- 
tionally invariant, 2.5D magnetohydrostationary configura- 
tion to the extended Grad-Shafranov equation. This section 
recapitulates how 2D equilibria balancing Lorentz force, 
gravity, and pressure gradients allow for axial (shear) flow, 
and are in principle realizable for three categories of solu- 
tions, where either density, entropy, or temperature take 
on a constant value on a single flux surface. We point out 
the essential role of straight field line coordinates, which 
must also be constructed from the obtained solution. In 
Sect. [3J we obtain the analytical equilibria valid for small 
gravity, where the outer flux surface is assumed to be cir- 
cular, and the internal, nested flux surfaces are then dis- 
placed circles, shifted downwards under the influence of 
gravit y. This is comp letely analogous to the Shafranov 
shift (Shafranov 1958) in the axisymmetric tokamak con- 
figuration induced by toroidal curvature. This can serve as 
a clear accuracy test for the numerically generated equi- 
libria. In Sect. HI we briefly explain the MHD equilibrium 
code FINESSE (jBelien et all 12002ft . This code is used to 
generate numerical equilibria. Finally, Sect. [5] then uses 
FINESSE and demonstrates how we obtain fourth-order ac- 
cu rate solut i ons by recovering the analytical case presented 
by IDungevi fl953) • It then continues to show the various 



types of solutions, how they obey the analytically predicted 
shift for small gravity, and then surveys those cases that are 
no longer tractable by analytic means. This includes a case 
where at sufficiently large gravity, multiply-layered promi- 
nence condensations can exist in a helical flux rope. We 
restrict these computations to cases for given outer cross- 
sectional shape resulting in only the internal solution. This 
restriction is needed because the follow-up analysis of di- 
agnosing the wave modes intrinsically depends on the pres- 
ence of closed, nested magnetic flux surfaces. The wave 
mod e analysis will be presente d in an accompanying pa- 
per (iBloklan d k Keppensll20lift . where analytic as well as 
numerical results on the MHD continua for the 2.5D equi- 
libria presented here are analyzed. 

2. Solar prominence equilibrium 

We begin our analysis by considering a translational 
symmetric solar prominence. For this kind of symmetry, 
Cartesian coordinates (x, y, z) are the natural choice. We 
note that for this symmetry, the equilibrium quantities de- 
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pcnd only on the poloidal coordinates x and y. The promi- 
nence equilibrium itself is modeled by the ideal MHD equa- 
tions 



<9v 

p— = -pv -Vv-Vp+jxB- pV$, 
at 



dp 

8t 
dB 

Th 

dp 

~dt 



— v • Vp — 7pV • v, 
V x (v x B) , 
-V • (pv) , 



(1) 
(2) 
(3) 
(4) 



where p, p, v, B, $, and 7 are the density, pressure, velocity, 
magnetic field, gravitational potential, and the ratio of the 
specific heats, respectively. The current density j = V x B 
and the equation V • B = has to be satisfied. We note 
that the presented ideal MHD equations are in the non- 
dimensional form. The equilibrium will be assumed to be 
time-independent. Furthermore, the thermal pressure and 
the density are related by the ideal gas law p = pT. From 
the equations V ■ B = and V • j = 0, it follows that the 
magnetic field and the current density can be written as 



B = e z xV(f + B z e z , 
j = -e 2 x V/ + j z e z , 



(5) 
(6) 



respectively, where, 2mp is the poloidal flux and the 
poloidal stream function I = B z . We restrict ourselves to 
prominence equilibria with a purely axial flow 



v = v z e. 



(7) 



In this case, Eqs. Q and (jl]) are trivially satisfied. The 
axial velocity v z is related to the electric field 



E 



v x B = v z Vip. 



(8) 



From the induction equation presented in Eq. ([3]) , it follows 
that v z — v z (ip). Since axial flows are routinely observed in 
prominence threads, we include the parameter v z in our 
analysis. Owing to the translational symmetry, it does not 
play any role in the ac tual force balance, but we s how in an 
accompanying paper (Blokland ~fc Keppensll20H how the 
velocity shear in the equilibrium modifies the eigenspec- 
trum. 

The momentum equation Eq. (TTJ) can be projected in 
three ways. The first one is in the axial direction, which 
shows that the poloidal stream function is a flux func- 
tion, i.e. / = I{ip). The second projection is parallel to 
the poloidal magnetic field resulting in two equations 



ip— const 



dp 
dx 

dp 

i/j=const 



d<S> 

d<5> 
~ P dy- 



(9) 



These two equations have to be satisfied simultaneously. We 
conclude from these that the pressure p = p{ip;x,y). The 
third and last projection is the one perpendicular to the 
poloidal magnetic field. This leads to the extended Grad- 
Shafranov equation 



^ dip dip 



(10) 



2.1. Further reduction for specific flux functions 

The two equations parallel to the poloidal magnetic field 
given in Eq. Q can be solved analytically under the extra 
assumption that either the temperature T, the density p, 
or the entropy S = pp 1 is a flux function. The assumption 
that the temperature is a flux function can be justified be- 
cause of the high thermal conductivity along the magnetic 
field lines. The resulting pressure can then be written as 

${x,y)~ 



p{ip;x,y) = p (ip)exp 



T(V>) 



(11) 



where po is the pressure for a static pure Grad-Shafranov 
equilibrium without gravity. The extended Grad-Shafranov 
equation in Eq. (fTU]l reduces to 



v 2 v = -I 



dl 

dip 



dpo 

dip 



Po$ dT 



T 2 dip 







exp 





(12) 



This extended Grad-Shafranov equation is a mino r gener- 
alization of the one presented in lPetrie et al.l ()2007t ). by al- 
lowing a fully general external gravitational potential. This 
case where temperature is a flux function is most relevant 
for the quiescent solar prominences, which are long-lived 
structures and hence will be temperature equilibrated on 
flux surfaces. The isothermal case forms a subset of this 
solution class. 

Another possibility, on MHD timcscales, is to assume 
that the density is a flux function. In this case, the pressure 
reads 

p(ip;x,y) =p Q {ip) [l 



where the quasi-temperature T p 



assumption, the extended 
Eq. (Tit)]) can be written as 



Po/p- 
Grad-Shafranov 



(13) 

Under this 
equation in 



V 2 ip 



dip 



dpo 

dip 



+ 



Po® dT p 
d^ 



7^2 
p 



$ 

1 



I 



(14) 



This case where density is a flux function will be shown 
in our companion paper to be susceptible to unsta- 
ble continuum modes, under q uantifiable stability crite- 
ria (|Blokland k. Keppens! 1201 it) . Hence, it might be rele- 
vant to the short-lived active region prominences, where 
these instabilities may relate to the sudden disappearance 
of filaments. 

The final option is to assu me that the entropy S is a flux 
function, as also discussed in lPetrie et al.l (|2007l ). We again 
generalize it slightly by adopting an arbitrary potential, for 
which the pressure reads 

7/(7-1) 



p(ip;x,y) =po(ip) 



1 - 



7-1 $(x,y) 



(15) 



7 T s {ip) 

and the extended Grad-Shafranov equation in Eq. (|10[) re 
duces to 



2/ T dI 



dpo 

dip 

1 - 



Po® dT s 

di/j 



rp2 
P 



7 - 1 $ 



7 T s 



(16) 



7 - 1 $ 
7 T s 



7/(7-1) 
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where the quasi-temperature T$ = Spq -1 . We note that 
all three cases can also be derived using the expressions for 
axisymmetric accretion t ori with a purely toroi dal flow, as 
presented in the paper bv lBlokland et al.f(|2007h . by setting 
R = Rq = 1 in their formulae. This case where entropy is 
a flux function is the natural one to consider when also 
plasma rotation would be incorporated. 

The equation for the density for all three cases can easily 
be derived by inserting the corresponding pressure equation 
into the equation for the momentum parallel to the poloidal 
magnetic field lines given in Eq. . The resulting equation 
is 



exp 
1 

1 - 



1 $ 



7 T s 



1/(7-1) 



for T = T(V>) 
for p = p(ijj) 

for S = S{t/j) 



(17) 

Here, the flux function p$ corresponds to the density of a 
related static equilibrium without gravity. 



2.2. Straight field line coordinates 

As we will adopt them for the actual stability analysis in 
our accompanying paper, we briefly discuss the 'straight 
field line' coordinates. These coordinates are an essential 
ingredient of an accurate stability analysis. For the conver- 
sion from the Cartesian (x, y, z) to straight field line coor- 
dinates (x 1 = ip,x 2 = d, x 3 = z), one needs the metric ten- 
sor and the Jacobian associated with the non-orthogonal 
coordinates in which the equilibrium field lines appear to 
be straight. Such a transformation is standard practice in 
MHD stability studies for laboratory tokamak plasmas. The 
metric elements </y and the Jacobian J are 

dr dr 

9lj = dx 1 ' dxl' (18) 



f =Vx l -Vx J , 
J = (V?/> x W ■ Vz) _1 , 



respectively. Here, the poloidal angle $ is constructed such 
that the magnetic field lines are straight in the z)-plane. 
The slope of these lines is a flux function 



dz 



ficldlinc 



B • Vz 
B • W 



= JI = qW), 



(19) 



where q is the safety factor. C omparing this expression with 
the one for tokamak plasmas (Wesson 20041) . one should re- 
alize that for tokamak plasmas the safety factor q is dimcn- 
sionless, while here t he factor q has a length dimension. As 
for tokamak plasmas (lGoedbloedlll975l : Ivan der Hoist et al.l 
2000; Blo kland et al.ll2007D . we introduce an expression for 
the poloidal curvature of the magnetic surfaces 

k p = -n • (t • Vt) = \ ( - 0^ ) JB® , (20) 
J \ 922/ 

where the unit vectors n = Vip/\Vip\ and t = H^/B^, 
and B$ is the poloidal magnetic field. The toroidal cur- 
vature Kt that is present in actual tokamak equilibria, is 
of course zero for a translational symmetric equilibrium. 
It is important to realize that the straight field coordi- 
nates can only be constructed when the solution ip(x, y) has 
been computed from the extended Grad-Shafranov equa- 
tion given in Eq. (ITU|) . 



3. Small gravity expansion 

In the previous section, we derived equations for the promi- 
nence equilibrium. In this section, we quantify the effect 
of gravity by means of a small gravity expansion. We will 
demonstrate in our companion paper that gaps will appear 
in the continuous MHD spectrum because of mode cou- 
pling, which is the result of the presence of gravity. This 
kind of expansion is similar to the small inverse aspect 
ratio e = a/i?o expansion for tokamak plasmas, where a 
and Rq are the minor radius of the pl asma and the ge - 
ometry axis of the tokamak, respectively (|Shafranov| [T958). 
Mathematically, both expansions, small gravity expansion 
and the small inverse aspect ratio expansion, are Taylor ex- 
pansions. In the remaining part of this paper, we assume the 
following gravitational potential in which the prominence 
is embedded 

= (x-x )g, (21) 

where xo is the location of the center of the last closed flux 
surface of the prominence and the gravity is represented by 
the constant g. 

The solar prominence equilibrium is expanded assum- 
ing that the gravity is small and that the outer flux surface 
is circular. Using these approximations, the fl ux surfaces 
can b e represented by slightly displaced circles (Shafranov 
(|1958f 0. which allows for the exploitation of non-orthogonal 
polar coordinates (r, 9, z), where r and 6 are the radius and 
the polar angle, respectively. Up to first order, we approx- 
imate 

%/}(x,y) = ip(r), 

x = x +rcos(6>) - A(r), (22) 
y = rsin(0), 



where A(r) is the Shafranov shift (|Shafranovl (|1958f l). which 
is expected to be in the downwards direction and caused by 
the gravity. As mentioned before, these polar coordinates 
are non-orthogonal and the associated metric elements are 



/i U w 1 - 2 A' cos(0), h 2 2 ~ r 2 
h 12 «rA'sin(60, ^33 = 1, 



(23) 



which means that the Jacobian J pa r[l — A'cos(#)]. 

Using the polar coordinates and expanding the extended 
Grad-Shafranov equation given in Eqs. (IT21) . (fT4l) . and (fl"6| 
up to first order leads to one equation for the equilibrium 
(0(1)) and one equation for the Shafranov shift (0(g)) for 
all three cases. The 0(1) equilibrium relation is 



dr 



Po 



B 2 e +B 2 z )] 



m 



0, 



(24) 



where Bg = 1// is the poloidal magnetic field expressed in 
polar coordinates, where the prime indicates the derivative 
with respect to the radius r. This is exactly as expected, 
since this equation merely expresses the force balance in a 
cylinder. From this equation, one can see that the pressure 
and the magnetic field components are of the same order 
O(l). This means that the plasma beta (3 = 2p/B 2 = 0(1). 
This is very different from tokamak physics where the pres- 
sure is two orders of magnitude lower than the toroidal mag- 
netic field component, which ensures that the plasma be ta 
there is of the order 0(e 2 ) (jvan der Hoist et all (|2000h ). 
For such tokamak plasmas, the Shafranov shift is out- 
wards due to the toroidic ity, pressure, and toroidal flow 
(|van der Hoist et al.l (|2000l )l 
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The equation for the Shafranov shift in prominence 
equilibria is 

which shows that the magnetic axis is shifted downwards 
because of the gravity. We note that for zero gravity or a 
constant density there will be no shift at all, again in accord 
with the cylindrical configuration then expected. 

As mentioned before, to perform the spectral formula- 
tion and compute the continuous MHD spectrum, the anal- 
ysis is done in straight field line coordinates. Exploiting the 
small gravity expansion, the straight field line metric ele- 
ments can be approximated by 



5ii 



1 

r 



[1 - 2A' cos(tf)] , 322 « r 2 [1 + 2A' cos(tf)] 



5 i2 « — (rA')'sin(tf), 033 = 1. 

(26) 

In deriving these expressions, we exploited that the safety 
factor q from Eq. (fT9"|) is a flux function and therefore the 
Jacobian J{= r/ip') has to be a flux function. This allows 
us to find a relationship between the polar angle 9 and the 
straight field line angle § 

0&-d + A'sin(i?). (27) 

4. The equilibrium code FINESSE 

4,1. Code basics: discretization 



The FINE S SE c ode was originally developed by 
iBelien et al.l (|2002l) and designed to solve the coupled 
generalized Grad- Shafranov equation toget h er with the al- 
;ebrai c Bernoulli equation (lHameiril (|1983l ); IZelaznv et al.l 
1993!): lGoedbloed et al.l (|2004h ). The design of the code is 
such that we could easily extend it by implementing the ex- 
tended Grad-Shafranov equations given in Eqs. (IT21 . (|T4")) . 
and <P~5|) . as realized in version 1.3 of the FINESSE code. 
These equations are solved for a given poloidal cross- 
section using the Picard iteration scheme in combination 
with a standard matrix solver. The boundary conditions 
are such that a fixed given boundary shape represents 
the last closed flux surface. As a discretization scheme, 
we exploited a finite element method in combination with 
the standard Galerkin method. As elements, we used 
isoparametric Hcrmitc elements. These elements ensure 
that the computed solution has the required accuracy 
needed for the stability analysis. 



4.2. Scaling for prominence equilibria 

The exploited scaling in FINESSE is such that all numerical 
quantities actually computed are quantified in 1) the profile 
variations of order unity, and 2) the amplitudes that deter- 
mine the relative strengths. The idea is that solving the gov- 
erning elliptic PDE numerically should exploit quantities of 
order unity as much as possible. This means for example 
that we ensure a scaling such that, e.g. the flux function ip 
varies from to 1. In the results section (Sect. [5}, we present 
either dimensionless quantities such as the local plasma /?, 
or quantify the pressure or the gravity parameter g in code 
units. Since we vary g in what follows, this may seem un- 
usual at first for solar filaments, because one can obviously 



not alter the solar gravitational field. We therefore here ex- 
plain the inherent scaling. 

For this purpose, some typi cal values for solar 
pro minence condit i ons t aken from lLabrosse et al.l {2010) 
and iMackav et alj ()2010f ) are as follows. Lengths are ex- 
pressed in units of the prominence-carrying magnetic loop 
radius. This length unit L can be varied from 100 km (if a 
single prominence thread is being modeled) to 10 Mm, the 
latter being an upper width quoted for solar prominence 
structures. Taking a typical temperature of 8000 K (a factor 
100 cooler than the corona) and an electron number density 
n e ~ 10 16 m~ 3 ,the density is then p ~ 1.673 x 10" 11 kg m~ 3 
and a reference pressure value is p ~ 0.00221 Nm -2 . We 
note that the sound speed can then be estimated to be 
15kms -1 . If we also use a typical value of magnetic field 
strength B ~ 10~ 3 T, the Alfven speed is 218kms~ 1 . 
For these parameters, the plasma beta can be as low as 
(3 = 0.0055, but we note that values of up to order 0.1 can 
be deduced from observations. The gravitational accelera- 
tion at a solar radius is of order g ~ 272 Nkg -1 , and a 
typical hydrostatic scale-height inferred from these values 
combined, yields A = p/pg ~ 485 km. When we quantify 
the parameter g in dimensionless units below, we actually 
mean the value of the dimensionless combination (with p,g 
the permeability) 



B 2 



0.0057. 



(28) 



We note that this dimensionless quantification of the 
strength of the gravitational field can meaningfully be al- 
tered by orders of magnitude, as it is directly proportional 
to the size of, and the density in the prominence, and in- 
versely proportional to the square of the field strength. 
Hence, larger values of the combination represent larger, 
denser prominences in weaker embedding magnetic fields, 
while the value given above is a reasonable reference value 
for a typical flux rope. Code units in essence use a dimen- 
sionalization based on L, p, and B, so that for instance the 
pressure is in units of B 2 /po ~ 0.795Nm~ 2 . This means 
that the reference pressure value quoted above is reached 
at p = 0.0028 in units of B 2 /p . 

5. Numerical results 

We now demonstrate that FINESSE can accurately com- 
pute prominence equilibria by comparing the n umerical 
solutio n with the analytical solution derived by iDuneevI 
(|l953l ). We then discuss two classes of equilibria. The 
first class represents cool prominences embedded in a hot 
medium and the second one in combination with proper 
chosen parameters represents multi-layered prominences. 



5.1. Accuracy test: the Dungey solution 

To demonstrate FINESSE accurately computes the promi- 
nence equilibria, we compare the numerical r esult with the 
analytical solution derived by IDuneevI (|1953h . The derived 
solution 



ip{x,y) = i> x 



ray 



1 



(29) 



where a = g/T, satisfies the extended Grad-Shafranov 
equation given in Eq. (|12j) with the extra assumption that 
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Fig. 1. The accuracy test of FINESSE with respect to the analytical solution derived bv iDungevl (|l953l) . The quantity 
L 2 and L°° are the average and maximum error between the Dungey solution and the numerical solution over the whole 
interior. 



the temperature T is constant. The solution itself is used 
to specify the last closed surface so FINESSE can compute 
the interior for various resolutions. The error between the 
numerical solution and the analytical one as a function of 
the number of grid points is shown in Fig.[TJ The quantities 
L 2 and L°° are the average error and the maximum error 
evaluated over the whole interior. The plot clearly demon- 
strates that FINESSE shows a fourth order convergence 
as it was designed to do. Furthermore, it also shows that 
moderate resolution is enough to obtain a highly accurate 
equilibrium solution. 

5.2. Cool prominence surrounded by a hot medium 

The first equilibrium class of solar prominences is the one 
where a prominence is embedded in a hot medium. These 
kinds of equilibria are specified by the following profiles: 



I 2 ^)=A(l-^\ 
p (V>) =AA 2 (1-^), 
T (VO =A 3 (i+9^ 2 -6^ 3 ), 
v z (if>)=0, 



(30) 
(31) 
(32) 
(33) 



where A^ = 0.1 and ^3 = 1. The amplitude A is the overall 
amplitude computed by FINESSE as part of the equilib- 
rium solution. The function To represents the (quasi) tem- 
perature of the three choices of the flux function, the den- 
sity, the entropy, or the temperature. The gravity g, scaled 
as explained earlier, has been varied from 0.001 up to 1.000. 
For most equilibria of this subsection, the temperature is 
assumed to be a flux function, except when stated differ- 
ently. Furthermore, we assume a circular boundary. The 
resolution adopted for the numerical computations are 101 



radial points and 129 points in the poloidal direction. The 
specifie d profiles for J 2 (t/>) and poW are of the same form 
use d bvlLow fc Hundhausenl (|l995fULow fc Zhand (|2004h . 
andlPgtrieeLalJ (120071) . The form of th e T a (j)) profile is the 
same as the one used bv iPetrie et al.l ([2007), which repre- 
sents a low temperature in the center while the temperature 
is high at the edge of the prominence. 

For the first equilibrium, we set the gravity to be g — 
0.001. The resulting pressure and plasma beta j3 = 2p/B 2 
are shown in Fig. [5] The plot shows clearly a high pressure 
in the center and a low pressure at the edge of the promi- 
nence. A similar observation can be made about the plasma 
beta, which is around 0.097 in the core and at the edge has 
the value of 0.014. The variation in the density, not shown, 
ranges from 0.006 at the edge up to 1.000 in the centre. 
The central temperature has a value of 0.048, while at the 
edge it is 0.780, which clearly resembles a cool prominence 
embedded in a hot medium. For the previously mentioned 
beta values, the magnetic field as well as the thermal pres- 
sure will play an important role in the stability analysis 
of our accompanying paper. In this stability analysis the 
safety factor q, shown in Fig. [21 is also an essential ingredi- 
ent. From fusion research, we know that around the q = 1 
and q = 3/2 surface, mode coupling may occur which in 
its turn can create gaps in the MHD continuous spectrum 
or even drive the continuous spectrum unstable. The plot 
of the safety factor shows that both surfaces exist in this 
cool prominence equilibrium. The same figure also shows 
the radial derivative of the Shafranov shift. This derivative 
was calculated in two different ways. The first applied the 
second of the equations given in Eq. (|2"2"|) and the ability 
of FINESSE to perform the transformation to straight field 
line coordinates. The second method numerically solved the 
Shafranov shift equation in Eq. (|25[) using Simpson's rule 
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Fig. 2. The two-dimensional pressure (red-scale) and plasma beta /3 = 2p/B 2 (contours) profile for a cool solar prominence 
surrounded by a hot medium with a gravity g = 0.001. The solar surface is below the figure. 
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Fig. 3. Left: the safety factor q at the mid-plane for a gravity g = 0.001. Right: the radial derivative of the Shafranov 
shift A(r) as a function of the radius r. 



(| Press et al.l (|1988|) ). The first method is valid for any value 
of the gravity parameter as long as the flux surfaces are cir- 
cles, while the second method can only be used for small 
gravity in combination with circular flux surfaces. The RHS 
plot of Fig. [3] shows excellent agreement between the two 
described methods. 

We then increase the gravity parameter by a factor of 
100, i.e. g — 0.100. The pressure and plasma beta are shown 



in Fig. SI Since the fact that the parameter g is 100 times 
stronger, the pressure maximum is shifted downwards as 
expected. The plasma beta has values similar to those in 
the previous case, its maximum and minimum values being 
0.105 and 0.014, respectively. The two-dimensional density 
profile, not shown here, has a similar behavior to that of 
the pressure profile. Its maximum is also shifted downwards 
as expected. Fig. [5] shows the safety factor q and the ra- 
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Fig. 4. The two-dimensional pressure (red-scale) and plasma beta (3 = 2p/B 2 (contours) profile for a cool solar prominence 
surrounded by a hot medium with a gravity g = 0.100. The solar surface is below the figure. 
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Fig. 5. Left: the safety factor q at the mid-plane for a gravity g — 0.100. Right: the radial derivative of the Shafranov 
shift A(r) as a function of the radius r. 



dial derivative of the Shafranov shift A(r). Comparing the 
safety factor with the one of the previous equilibrium, one 
notices that there is hardly any difference. The two methods 
for computing the Shafranov shift now show a small differ- 
ence. The second method, which using the Shafranov shift 
equation in Eq. (|25|) . underestimates the radial derivative 
for all radii. 



The last cool solar prominence we discuss is the one 
for which the dimensionless gravity parameter is g — 1.000. 
The gravity is 1000 times stronger than in the first discussed 
equilibrium and 10 times stronger than the last one. The 
downward shift of the pressure is even stronger as shown in 
Fig. [6] For this strong g case, the plasma beta varies from 
0.013 at the edge up to 1.107 at the core of the plasma. 
The latter value is larger than one would expect from ob- 
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Fig. 6. The two-dimensional pressure (red-scale) and plasma beta (3 = 2p/B 2 (contours) profile for a cool solar prominence 
surrounded by a hot medium with a gravity g = 1 .000. The solar surface is below the figure. 
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Fig. 7. Left: the safety factor q at the mid-plane for a gravity g = 0.500. Right: the radial derivative of the Shafranov 
shift A(r) as a function of the radius r. 



servations. However, it is straightforward to compute an 
even more realistic equilibrium by adjusting the coefficient 
Ai of the pressure equation in Eq. (|31j) . For this equilib- 
rium, the safety factor and the Shafranov shift are shown 
in Fig. [7] Owing to the strong gravity the safety factor is 
very different from its value for the two previous equilib- 
ria, particularly in terms of the multiple q = 1 surfaces. In 
addition, the plot of the Shafranov shift shows a large dis- 



crepancy between the two methods, because of the strong 
gravity. The first method of determining the Shafranov shift 
is superior in terms of accuracy and numerical cost. 

To illustrate the difference between the chosen flux func- 
tion, we also computed the cases where the density or the 
entropy is a flux function for gravity g — 1.000. For these 
choices, the pressure and plasma beta are plotted in Fig. [SI 
When density is a flux function, the plasma beta varies 
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Fig. 8. The two-dimensional pressure (red-scale) and plasma beta /3 = 2p/B 2 (contours) profile for a cool solar prominence 
surrounded by a hot medium with a gravity g — 1 .000 for the case that the density (left) or the entropy (right) is a flux 
function. The solar surface is below the figure. 



from O.Off up to 0.254, and for the entropy case the range 
is O.Of 2 up to 0.395. Both plots are clearly distinct from the 
case where the temperature is chosen to be a flux function. 
From the equilibrium viewpoint, all three are realizable, but 
high resolution observations will be required to help us to 
select the proper flux function for a particular filament. 

5.3. Multi-layered prominences 

The second class of equi li bria, based on the description 
presented by iPetrie et al.l (|2007l) , allows for multi- layered 
prominences. This class is specified by the following ex- 
pressions 

I 2 (^)=A(1-^), 
po(i/>) = AA 2 (l - tp - 18i/j 2 + 56±V 3 



(34) 
(35) 



59V 4 + 20fiV 5 ), 



ToW=A 3 , 

u.W0 =0, 



(36) 
(37) 



where A 2 — 0.1 and A3 = 1. In addition, the overall ampli- 
tude A is computed by FINESSE. The gravity g is varied 
from 0.001 up to 5.000. As for the previous class of cool 
prominences, we assume that the boundary is circular and 
that the temperature is a flux function. The function To 
has the same meaning as in the previous subsection, and 
represents the (quasi) temperature of the chosen flux func- 
tion. As for the previous class, the numerical computations 
are done for 101 and 129 points in the radial and poloidal 
direction, respectively. 

Starting with small gravity g — 0.001, Fig. [S] shows the 
two-dimensional pressure and the plasma beta j3 = 2p/ B 2 . 



The latter quantity varies from 0.001 at the edge to 0.097 
at the center of the prominence. Once more, we empha- 
size that for these plasma beta values the magnetic field as 
well as the pressure play an important role in the stability 
analysis of our accompanying paper. The pressure shows a 
ring structure. By increasing the parameter g, we show be- 
low that one can create a double-layered prominence. The 
safety factor q and the radial derivative of the Shafranov 
shift are plotted in Fig. [TU] For small g, the safety factor 
q has two flat regions that correspond to the ring-like cav- 
ity of low pressure. Furthermore, note that the q = 1 and 
g = 3/2 surfaces both exist for this equilibrium. As for cool 
prominences, this will mean that gaps or instabilities may 
occur in the MHD continuous spectrum. The figure also 
shows that the two methods for computing the Shafranov 
shift, as discussed in the previous subsection, show excellent 
agreement. 

For the next equilibrium, we have set the gravity param- 
eter g = 1.000, which is 1000 times stronger than in the case 
above. This example shows the onset of the double-layered 
pressure profile as can be seen in Fig. 1111 If one closely 
looks at the plot, one still sees the ring structure but be- 
cause of the stronger gravity compared to the small grav- 
ity 9 = 0.001 case this structure has almost disappeared. 
The plasma beta ranges from 0.00069 in the low pressure 
regions up to 0.117 in the high pressure regions. For this 
equilibrium, we also computed the safety factor q and the 
Shafranov shift. Both quantities are presented in Fig.[T21 A 
comparison between the safety factor for this equilibrium 
and that for the small gravity g = 0.001 case finds hardly 
any differences. The only real difference is the maximum of 
the q profile, which is shifted slightly downwards because of 
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Fig. 9. The two-dimensional pressure (red-scale) and plasma beta (3 = 2p/B 2 (contours) profile for a double-layered solar 
prominence with a gravity g = 0.001. The solar surface is below the figure. 
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Fig. 10. Left: the safety factor q at the mid-plane for a gravity g = 0.001. Right: the radial derivative of the Shafranov 
shift A(r) as a function of the radius r. 



the stronger gravity. The plot of the Shafranov shift shows 
that a small difference between both means of determining 
the shift is clearly evident. The second way of applying the 
Shafranov equation given in Eq. ([251 again underestimates 
the radial derivative. 

The increase in the gravity parameter to g — 5.000 cre- 
ates a clear double-layered pressure structure as can be seen 
in Fig. [T3] A similar structure is found for the density. The 



plasma beta varies from 0.000031 up to 1.994. The latter 
value is higher than we expect based on observations. As 
mentioned before, it is straightforward to compute a more 
realistic value by changing the coefficient of the pro- 
file pressure equation in Eq. (|35p . The safety factor of this 
double-layered prominence, shown in Fig. 1141 varies widely. 
It no longer contains a q = 3/2 surface as in the previ- 
ous two cases of this equilibrium class. Furthermore, we 
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Fig. 11. The two-dimensional pressure (red-scale) and plasma beta /3 — 2p/B 2 (contours) profile for a double-layered 
solar prominence with a gravity g — 1.000. The solar surface is below the figure. 
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Fig. 12. Left: the safety factor q at the mid-plane for a gravity g = 1.000. Right: the radial derivative of the Shafranov 
shift A(r) as a function of the radius r. 



note that near the magnetic axis (around x = —0.5) it the varying twist of the magnetic field lines. In our accom- 
almost creates new q = 1 surfaces. The magnetic axis is panying paper, we will analyze the way in which this more 
the location where the poloidal magnetic field is zero. The realistic field configuration modifies the MHD eigenspec- 
Shafranov shift is also shown in Fig. [TU As expected, the trum of the flux-surface localized continuum modes, 
two methods for determining the shift no longer agree. As in 
the previously discussed case, the Shafranov shift equation 
in Eq. (pj)) underestimates the radial derivative. Figure [T51 
presents this prominence configuration in 3D, illustrating 
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Fig. 13. The two-dimensional pressure (red-scale) and plasma beta (3 — 2p/B 2 (contours) profile for a double-layered 
solar prominence with a gravity g = 5.000. The solar surface is below the figure. 
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Fig. 14. Left: the safety factor q at the mid-plane for a gravity g = 5.000. Right: the radial derivative of the Shafranov 
shift A(r) as a function of the radius r. 



6. Conclusions 

We have presented the equations describing translational 
symmetric solar prominence equilibria. These equations 
have been specialized for three choices of the freedom in 
flux function in which the temperature, density, or entropy 
can be flux-dependent. For small gravity and circular cross- 
section, we have expanded these equations to derive the 



equation for the equilibrium of a cylindrical plasma and 
the equation for the Shafranov shift. 

The MHD equilibrium code FINESSE has been used 
to compute accurate prominence equilibria. A compari- 
son between the numerically generated equilibrium and the 
Dungey solution shows fourth order convergence, and mod- 
erate resolution is needed to obtain an accurate solution. 
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Fig. 15. A three-dimensional impression of the double-layered solar prominence for gravity parameter g — 5.000. Shown 
are an isosurface of density, at a value showing the layering of the condensation, with the cross-sectional view shown 
midway this 3D impression. The magnetic field structure is visualized by drawing selected field-lines, demonstrating the 
varying twist of the helical configuration. The held lines are colored by the magnetic field strength. The solar surface is 
below the figure. 



We have considered two equilibrium classes, cool promi- 
nences surrounded by a hot medium and double-layered 
prominences, for different values of the dimensionless grav- 
ity parameter (in a range corresponding to solar filament 
cases) with plasma beta of the order 0.1. For both classes, a 
comparison between the numerical solution and the derived 
equation for the Shafranov shift has been made. This com- 
parison shows excellent agreement for small gravity, while 
for strong gravity only method 1 is able to quantify the 
shift. 

For cool prominences, the results show that by increas- 
ing the gravity, the location of the maximum pressure shifts 
in the downwards direction. This downward shift is largest 
when the temperature is a flux function. Furthermore, in a 
strong gravitational potential the three choices of the cho- 
sen flux function show large deviations from each other in, 
for example, the pressure and density. This means that ob- 
servations could indicate what the most suitable choice of 
the flux function should be. 

The results of the double-layered prominences have re- 
vealed that a double-layered structure in pressure and den- 
sity can be created in actual filament configurations. The 
relative strength of the gravitational potential must be 
sufficient: if the potential is weak, a ring structure ap- 
pears, where a cavity surrounds the prominence. As for 
cool prominences, the location of the maximum pressure 



and density is shifted increasingly downwards if the gravity 
importance is increased. 

In. our accompanying paper (Bloklan d fc Keppensl 

1201 ll ). the stability properties of these equilibria will be an- 
alyzed, with special attention to the continuous MHD spec- 
trum. Owing the presence of gravity, gaps or even instabili- 
ties may appear in this continuous spectrum. Furthermore, 
inside these gaps new global modes may occur, which pro- 
vide us with important information about the internal 
structure of the prominence. Before investigating the pos- 
sible appearance of global modes, a detailed analysis of the 
continuous spectrum will be required. 
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